
#ifndef __wrapper_h__
#define __wrapper_h__

#include "petsc_vector.h"
#include "petscsnes.h"

  // Forward declaration
class NeutronicsSystem ;

namespace Wrapper
{
  typedef PetscErrorCode(*orthonorm_function)(Vec, Vec) ;

  void initialize_wrapper (const PetscVector<Real>& ) ;

  void clear_wrapper () ;

    // utility function to copy a smaller vector to a larger vector and vice-versa
    // useful when solving simple eigen-problem by Newton
  PetscErrorCode petsc_copy_vector(Vec xin, Vec xout, bool is_input_smaller=false) ;

    // y = inv(L)*P*x; Mat represents the operator = inv(L)*P
  PetscErrorCode simple_eigen_operator(Mat, Vec, Vec) ;

    // y = lambda*x - inv(L)*P*x
  PetscErrorCode simple_eigen_residual(NeutronicsSystem *, Vec, Vec ) ;

    // y = L*lambda*x - P*x
  PetscErrorCode generalized_eigen_residual(NeutronicsSystem *, Vec , Vec ) ;

    // Preconditioner for the simple eigen oeprator: useful for solving with ST = SIPI
  PetscErrorCode simple_eigen_preconditioner_solve_operator(PC pc, Vec x, Vec y) ;

    // y = L*lambda*x - P*x
  PetscErrorCode generalized_eigen_residual(NeutronicsSystem *, Vec , Vec, double ) ;

    // y = L*x; If Mat is shell, return operator application; If Mat is valid, return product
  PetscErrorCode multiply_loss_matrix(Mat , Vec , Vec ) ;

    // y = P*x; If Mat is shell, return operator application; If Mat is valid, return product
  PetscErrorCode multiply_gain_matrix(Mat , Vec , Vec ) ;

    // Generic wrapper that finds y = J*x for the simple/generalized eigenvalue problem
    // This assumes that Keff() has been updated by residual or some other call explicitly
  PetscErrorCode multiply_jacobian_matrix(Mat, Vec, Vec) ;

    // Generic wrapper that finds y = J*x for the simple/generalized eigenvalue problem
    // This calculates Keff() from the solution x given and then finds J*x
  PetscErrorCode multiply_true_jacobian_matrix(Mat, Vec, Vec) ;

    // y = J*x for the simple eigenvalue problem
  PetscErrorCode multiply_simple_jacobian_matrix(NeutronicsSystem *, Vec, Vec) ;

    // y = J*x for the generalized eigenvalue problem
  PetscErrorCode multiply_generalized_jacobian_matrix(NeutronicsSystem *, Vec, Vec) ;

    // If possible, implement a function to evaluate the diagonal of Shell matrix
  PetscErrorCode shell_matrix_get_diagonal(Mat, Vec);

    // If possible, apply a shift to the main diagonal of the operator Shell matrix
    // This is important for shift-invert strategies to work correctly.
  PetscErrorCode shell_matrix_apply_shift(Mat, PetscScalar);

    // Based on the wrapper routines defined for the shell matrices in the context of
    // eigenvalue computation, copy the Shell matrix by making sure all properties
    // of given matrix are propagated; Very similar to Duplicate routine since
    // there are no entries stored for Shell matrix.
  PetscErrorCode shell_matrix_copy (Mat , Mat , MatStructure) ;

    // Based on the wrapper routines defined for the shell matrices in the context of
    // eigenvalue computation, duplicate a Shell matrix by making sure all properties
    // of given matrix are propagated.
  PetscErrorCode shell_matrix_duplicate (Mat , MatDuplicateOption , Mat *) ;

    // Define a dummy method to set Matrix entried to zero
    // Since most matrices are Shell type, nothing much to do here !
  PetscErrorCode shell_matrix_zero_entries (Mat a) ;

  PetscErrorCode precondition_loss_matrix_solve(PC , Vec , Vec ) ;

    // Solves the Loss matrix through block-gauss-seidel method ; 1 sweep
    // serves as a preconditioner to Loss matrix inversion
  PetscErrorCode precondition_loss_matrix(NeutronicsSystem *, Vec , Vec ) ;

    // Generic wrapper that multiplies Preconditioner matrix representing the Jacobian for
    // simple/generalized eigen value problem
  PetscErrorCode preconditioner_multiply(Mat , Vec , Vec ) ;

    // Generic wrapper that solves Preconditioner matrix representing the Jacobian for
    // simple/generalized eigen value problem
  PetscErrorCode preconditioner_solve(PC, Vec , Vec ) ;

    // Solves the Preconditioner matrix that represents the Jacobian for
    // simple eigen value problem
  PetscErrorCode simple_preconditioner_solve(NeutronicsSystem *, Vec , Vec ) ;

    // Solves the Preconditioner matrix that represents the Jacobian for
    // generalized eigen value problem
  PetscErrorCode generalized_preconditioner_solve(NeutronicsSystem *, Vec , Vec ) ;

  void set_orthonormalize_reference(orthonorm_function ) ;

}

extern "C" PetscErrorCode eigen_residual (SNES snes, Vec x, Vec r, void *ctx) ;

extern "C" PetscErrorCode eigen_jacobian (SNES , Vec x, Mat *jac, Mat *prec, MatStructure *msflag, void *) ;

extern "C" PetscErrorCode snes_monitor(SNES , PetscInt , PetscReal , void *) ;

#endif // #define __wrapper_h__
